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Time-evolving bubbles in two-dimensional Stokes flow 


Saleh Tanveer 1 and Giovani L. Vasconcelos 
Department of Mathematics 
Ohio State University 
Columbus, OH 43210 


Abstract 

A general class of exact solutions is presented for a time evolving bubble in a two- 
dimensional slow viscous flow in the presence of surface tension. These solutions can describe 
a bubble in a linear shear flow as well as an expanding or contracting bubble in an otherwise 
quiescent flow. In the case of expanding bubbles, the solutions have a simple behavior in the 
sense that for essentially arbitrary initial shapes the bubble will asymptote an expanding 
circle. Contracting bubbles, on the other hand, can develop narrow structures (‘near-cusps’) 
on the interface and may undergo ‘break up’ before all the bubble-fluid is completely re- 
moved. The mathematical structure underlying the existence of these exact solutions is also 
investigated. 


Research was supported in part by the National Aeronautics and Space Administration under NASA 
Contract No. NAS1-19480 while the author was in residence at the Institute for Computer Applications in 
Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23681. 
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1. Introduction 


The study of the deformation and breakup of drops and bubbles in a slow viscous flow 
is of practical significance to many physical processes such as the rheology of emulsions and 
the mixing in multiphase viscous systems. Following the pioneering work by G. I. Taylor 
(1932, 1934), there has been a great deal of both theoretical and experimental research on 
the subject. The reviews by Acrivos (1983) and Rallison (1984) summarize the state of 
affairs in the early eighties. In many of these early studies the term 'drop breakup’ usually 
does not refer to the fragmentation of a drop, but to the non-existence of a steady solution 
when the applied shear strength exceeds some critical value. In the past decade, however, 
there have been a number of mainly experimental and computational investigations of the 
actual dynamics leading to breakup. These have recently been reviewed by Stone (1994). 

Analytical solutions for the time evolution of a general three-dimensional drop or bubble 
in a slow viscous flow do not appear amenable to currently known techniques. The simplified 
case of two-dimensional bubble flows, on the other hand, is analytically tractable through 
complex variable methods; and their study might shed some light on important qualitative 
aspects of axisymmetric 3D flows. In this vein, Richarson (1968) has obtained exact solutions 
for an inviscid 2D bubble in a linear shear flow, while Buckmaster & Flaherty (1973) has 
found approximate solutions for a drop with the same viscosity as the ambient fluid. 

There have also been recent investigations of 2D Stokes flows with time-dependent free 
boundaries. In a series of papers, Hopper (1990, 1991, 1992, 1993) has found several analytic 
solutions for the motion of a blob of viscous fluid driven by surface tension. Among these 
are solutions describing the coalescence of two cylinders (Hopper 1990) and the coalescence 
of a cylinder with a half-plane (Hopper 1992) — two problems of interest in viscous sinter- 
ing. Richardson (1992) has recently reviewed the mathematical structure of these solutions. 
Howison and Richardson (1993) have also obtained exact solutions for the case of a 2D 
viscous blob with n-th fold symmetry (i.e., invariance under rotation by 2x/n , for n > 1) 
in the presence of suction. In the absence of surface tension their solutions develop cusp 
singularities on the interface (and hence cease to be physically meaningful) before the fluid 
is completely removed. A nonzero surface tension, however, allows for total removal of fluid. 

The main aim of the present paper is to report a general class of exact solutions for a 
time-evolving bubble in a 2D Stokes flow. Our solutions include, for instance, the motion of 
a bubble of essentially arbitrary initial shape in a linear flow or an expanding/contracting 
bubble in an otherwise quiescent flow. In the former case, one is able to follow the transient 
motion of the bubble as it approaches the steady state. In an earlier paper (Tanveer k 
Vasconpelos 1994a), we briefly discussed the case of specific solutions that show a contracting 
bubble developing topological singularities (‘breakup’) before the bubble shrinks to zero 
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size. Other solutions show that cusp can form in the presence of certain symmetries when 
surface tension is neglected-in the presence of surface tension, however, cusp formation is 
inhibited until all the bubble-fluid (‘air’) is extracted. It was also pointed out that contrary 
to the expectations from a Hele-Shaw flow analogy recently made in the literature (Howison 
and Richardson 1993), a circular bubble expanding in a 2D Stokes flow is stable, while a 
contracting one is not. Here in this paper, we report a more complete set of such solutions 
for expanding/contracting bubbles, present detailed analysis of the ‘regularization’ effects 
provided by surface tension, and discuss the mathematical structure that, among other 
properties, guarantees the existence of exact solutions for a rather general class of initial 

conditions. 

The paper is organized as follows: in §2 the problem is mathematically formulated in 
terms of a conformal mapping from the interior of a unit circle to the flow domain. It should 
be noted here that the formulation of the problem of a bubble in a 2D Stokes flow parallels 
that of a 2D viscous drop, as presented by Hopper (1990) and Richardson (1992). In fact, 
a shorter presentation in §2 could have been achieved by quoting some of their formulas. 
However, since the two problems differ in a number of aspects and since there are differences 
in the derivations, we felt that a complete formulation was warranted here. In §3 we discuss 
certain global properties of the conformal mapping that underlie the existence of exact 
solutions for the problem. Many of the arguments presented in this discussion transcends 
the specific details of the problem in question and have been found to apply to other 2D 
free-boundary problems as well (Tanveer 1993). Accordingly, this section is likely to appeal 
to researchers interested in finding exact solutions to free-boundary problems that can be 
conveniently cast in terms of conformal mapping. For those more interested m the concrete 
results for a 2D Stokes bubble, however, it can be skipped without any loss of continuity. In 
§4 a general class of exact solutions of polynomial type is presented, while the following two 
sections give details for two specific cases. First, in §5 the problem of a bubble placed in a 
linear flow is considered for two particular flow arrangements of relevance to experiments: 
§5.1 describes a bubble in a pure straining flow, whereas §5.2 focuses on a simple shear flow. 
The case of expanding/contracting bubbles in an otherwise quiescent flow is then addressed 
in §6. Our conclusions and main results are summarized in §7. 

2. Mathematical Formulation 

We consider the problem of a bubble placed in a two-dimensional slow viscous flow. The 
fluid inside the bubble has a negligible viscosity and is at a constant pressure, which is 
chosen to be zero without loss of generality. The fluid outside the bubble has a viscosity 

and is incompressible. Under the assumption of no inertial effects, gravitational or other 
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body forces, the fluid motion is governed by the Stokes equation and the incompressibility 
condition: 


fiV 2 u = -Vp (1) 

Vu = 0, (2) 

where u(.r,p) is the fluid velocity and p the pressure. 

On the bubble boundary we must ensure continuity of the shear stress and satisfy the 
requirement that the jump in the normal stress across the interface equals the product of the 
surface tension a and the curvature k. These two stress boundary conditions can be written 
as 


-prij + 2p e jk n k = ctk rij, ( 3 ) 

where the indexes j and k take on values 1 and 2 (the Einstein summation convention is 
used in the above), and n\ and n 2 are the x and y component of the outward normal unit 
vector. Here e jk are the components of the rate of strain tensor and are given by 


e ifc 


1 / duj du k 

2 V&rjt ^ dx 


-)■ 

3 7 


( 4 ) 


In addition, we must also satisfy the usual kinematic condition that the normal velocity V„ 
of a point on the bubble surface is equal to the normal component of fluid velocity at that 
point, that is, 


u n = V n . (5) 

To completely specify the problem we need to prescribe appropriate boundary conditions 
at infinity. We suppose in general that the bubble is placed in a linear flow and its area is 
changing at a prescribed rate m, which in general can be a time dependent function. More 
specifically, we require that far away from the bubble the fluid velocity behaves as 

jjl x 

u ~ T-x + + 0(i/|x:| 2 ), for |x| -► oo, (6) 

where 


r = I ( °° Po-Uo\ , . 

2 \ A + wq — J ' ^ 

Here u> 0 is the vorticity of the external flow, while a 0 and /3 0 characterize its strain. 

It should also be noted at this stage that the problem stated above can be recast in terms 
of nondimensional quantities if we rescale velocities by a/p, pressure by a/R , and the length 
and time scales by R and Rp/cr , where ttR 2 is the bubble initial area. In these units, the 
dimensionless parameters characterizing the problem are 

r = T and m' = m. (8) 

<7 (7 R K 1 
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We prefer however to use dimensionful quantities throughout the rest of the paper, except 
where noted otherwise. This is because a different nondimensionalization is appropriate 

when <7 = 0. 

As is well known (Lamb 1945), the problem of two-dimensional Stokes flow can be con- 
veniently formulated in terms of a stream function ip{x,y), defined as 


dip 

u i = -5-, u 2 

dy 


dip 

dx 1 


( 9 ) 


so that ip relates to the fluid vorticity w through 

V 2 0 = -u ( 10 ) 


and obeys the biharmonic equation 

VV = 0. 


( 11 ) 


Alternatively, one can formulate this problem in terms of a stress function <£(x,y) (Muskhe- 
lishvili 1963), defined via 

( 12 ) 


V 2 ^ = 

a 


so that also obeys the biharmonic equation. Here the time dependence of 0 and <f> has 
been omitted for notational convenience. 

Next we introduce the quantity W(z, z) = <p{x, y)+iip(x, y), where z = x + iy and the bar 
denotes complex conjugation. Then according to the Goursat representation for biharmonic 
functions (Carrier, Krook, and Pearson 1966), W{z,z) can be written as 


W{ Zl z) = zf{z)+9(z)i 


( 13 ) 


where f(z) and g(z) are analytic functions in the fluid region. All the physically relevant 
quantities can now be expressed in terms of the functions f(z) and g{z). After a little 
algebra, one can easily establish the following identities: 


p 

ZUJ 

= 4 f{z) 

(14) 

P 

(15) 

U\ + IU2 

= — fi z ) + z f'(z) + d'(z), 

e n + 

= zf"{ z ) + g"(z). 

(16) 


where prime indicates derivative and / denotes the conjugate function. f(z) f(z) (and 
similarly for g). The functions f(z) and g\z) must also satisfy appropriate boundary condi- 
tions, as described below. 

First consider the boundary conditions at infinity. From (14) it follows that 


/(*) 


Po o(0 


-!W 0 


z + C{t) + 0(l/z), as |z|-»oo, 


(17) 
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where the functions Poo(t ) [the pressure at oo] and C{t) are to be determined later. Now 
using (6), (15), and (17), one finds 

d '( z ) ~ ^ ( a ° ~ *A>) ^ + C(t) + + 0(1/ ^ 2 ), as |z| — > oo. (18) 

Note that according to (14) and (15) the choice of C(t ) does not affect the velocity and 
pressure fields a specific choice however will be made below for convenience. 

If we now define 

N ~n j + in 2 = i(x s + iy») = iz s = ie ,e , (19) 

where s is the arclength traversed in the clockwise direction and 0 is the angle between the 
tangent and the real positive x axis, then the two stress conditions given in (3) can be written 
as one complex equation: 


— pN + 2p(eu + ie u )N = crnN. 


( 20 ) 


Using (14), (16), (19), and the fact that k = —9 S , we find after a straightforward calculation 
that (20) is equivalent to 


dH(z,z) , = dH{z,z) 

Z s O "f" z s o- 
OZ Oz 


( 21 ) 


where r = a j p and 

H{z,z) = f(z) + zf'(z) + g'(z). ( 22 ) 

Equation (21) can then be integrated once, yielding the following condition on the bubble 
surface: 

f ( z ) + z f \ z ) + g{z) = ~^2 Z s- (23) 

Here, without loss of generality, the constant of integration has been set to zero. [This 
corresponds to a specific choice of C(t).] From (15) and (23), it also follows that on the 
bubble surface: 

ui + iu 2 = -i^z s - 2f(z). (24) 

Next, we consider the conformal mapping z((, t) that maps the interior of the unit circle 
in the C plane to the fluid region (i.e., the exterior of the bubble) in the z plane, such that 
the = 0 corresponds to the point z — co. We thus write 


*(C,0 = ^ + A(C,<)» (25) 

where a{t) can be chosen real and negative in view of the additional freedom of the Riemann 
mapping theorem. Here h((,t ) is assumed analytic in |(| < 1 and such that / 0 in |([ < 1 
at least for some period of time. 
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Now recall that the kinematic condition (5) can be written in complex notation as 

Im {[ Zt — («i + * 1 x 2 )] = 0, (26) 

where Z{u,t) = z{e ip ,t). Note also that 

(27) 


_ iCfc 
z <l 


z. = T— t , on Id = 1. 


Inserting (24) and (27) into (26) and using the fact that Z v = i(z ( on ( = e", we obtain the 
following condition on |d = 1: 


Re 


z t + 2 


C z < 


2|*d’ 


(28) 


where r = crj n and 

F((,t) = f(z((,t),t). (29) 

For later use, we note that that in view of (25) the asymptotic condition (17) implies that 

F{(,t) ~ tksdtz + C(t) + 0(0, as |z| — »■ 0. (30) 

We then notice that the quantity within square brackets in (28) is an analytic function of 
C for ICI < 1. [Note that the simple poles at ( = 0 in both the numerator and denominator 
cancel out.] From Poisson formula (Carrier, Krook, and Pearson 1966), it follows that for 

ICI < 1: 

%t + 2F(C 

where A is a real constant and 

1 (^ = 7 - 1 


= d^(Cd) + iF] Z{, 

(31) 

dc r c + cl 

1 

(32) 

vTN j 

f 

MC',*)I' 

► 0, we readily obtain 


1 

— L*>0 , 


(33) 

2 

-2n f/(0,0 + 

=], 

(34) 


where the dot denotes time derivative. Equation (31) can thus be thought of as an evolution 
equation for the mapping function z{(,t). Indeed, in the next section we shall use this 
equation to derive several general properties concerning the motion of the singularities of 
2 (C,0 For practical purposes, however, it is more convenient to view (31) as determining 
F((,t) [and hence the pressure and vorticity fields] in terms of the conformal mapping. 
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Now if we define 


G (f.O =$'(*(C,<)A (35) 

then the function (?(£,<) can also be expressed in terms of z((, t) as follows. Inserting (27) 
into (23), taking the complex conjugate, and using the fact that C = 1/C on \(\ = 1, we then 
find 

‘V.o 


G(C. () = -F(C\ 0 - K <-\ J 


2 <(C,<) ' <36) 

which is originally valid on |C| = 1 and therefore elsewhere by analytic continuation. Using 
(31) to eliminate from the equation above, we obtain after some simplification 


G«,t) = ^(C 1 


+ 



Ch(C,t)~ 

1 


1 + 


C 2 «(C? 0 


2 <(C,<) 


[/(C,i) + *A'] 


[I(( y t) + zK} + -z t (C-\t). 


[In deriving this, we have made use of the relation 

+ KC\t) = 


(C.^UlC-V) 


(37) 


(38) 


which follows from the fact that Re /((,£) = r/2|z<;(C, t)| on |C| = 1.] 

We point out at this stage that (37) not only gives G((,t) in terms of z((,t) but it also 
determines the evolution of the mapping function itself. In fact, the requirement that the 
right hand side of (37) is analytic in |C| < 1, except for a known pole at C = 0, can be used to 
determine the time evolution of the parameters in our solutions, as we will see in §4. Before 
doing that, however, we would like to describe in the next section the general mathematical 
properties that underlie the existence of a rather broad class of exact solutions. 


3. General Properties of Singularities 


In this section we shall derive several global properties of the mapping function z(( , t ) for 
arbitrary initial data £(C,0) analytic in |C| < 1 (except for a simple pole at C = 0) and with 
2 c(C>0) 7^ 0 in |C| < 1 as well. In particular, we shall be interested in the behavior of the 
singularities of in |(,| > 1. We begin by analytically continuing (31) into the exterior 

of the C unit circle. Through a standard procedure of contour deformation, one finds that 
for |£| > 1: 


*. + 2F(C, () = <[/«, i) + ,K}z ( + r j fr ) C ’° . 

z c, (C 

Now on taking the complex conjugate of (36) on the unit ( circle, it follows that 

<j(c\f) = -F«,o - ,«,.) M^ + i # ( c.o . 

2 c(C .*) 2 ^ /2 (C-',0 


(39) 


(40) 
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on | C | = 1 and elsewhere by analytic continuation. In view 

of (40), it is clear that (39) can 

be written as 

Zt = qiz^ + q 3 z + < 72 , 

(41) 

where 

?.«,() = C [<«.*) + **"1 

(42) 


,, ,, J<Sfr'A 

(43) 



(44) 


Note from the definitions (32) and (42)-(44) that as long as a solution z{(,t) exists for 
which z c (C, t) is nonzero and analytic in |(| < 1 each of ( -1 <7i, ?2 and <73 will remain analytic 
in |C| > 1. This implies that the corresponding Laurent series on |(| = 1 will contain only 
nonpositive powers of (. For instance, one has that on the unit circle 

0= E «».(0 < 45) 

n=— oo 

Similar Fourier series representations exist for and e~ ,l, qi(e xy , t). 

In order to further elucidate the properties of the evolution equation (41) in |d > 1, it is 
convenient to introduce the following projection operator Vn acting on the class of functions 
v(v) with a convergent Fourier representation 

OO 

»(*)= E (46) 

n= — co 

where Vn for N > 0 is defined by 

OO 

IP»»1M= E (47) 

n=N+l 

On applying the operator Vn to (41), it follows that the projection = Vn z(e x \t) 

satisfies the following evolution 

H t = V N \-ie~ iv q\(e iv ,t)Hj\ + V N [q 3 {e w ,t)H] . (48) 

It follows from (48) that if H( i/,0) = 0, then H(v,t) = 0 for all times. This fact implies in 
turn that if z(C,0) has no singularities in |(| > 1, except for a “pole” of order N at infinity, 
i.e., z(C,0) ~ Constant • ( N as ( — *■ oo, then z((,t) will have the same property for t > 0. 
Since z{(,t) is a conformal map with only one possible singularity in |£| < 1 at £ = 0, it 
follows that if we start with an initial condition of the form 

2(C,o) = ^ + MC»). (49) 

8 


i 1 1 


with /i((,,0) a N-th order polynomial, then for later time we will have 

*«,*) = + (50) 

where h((, t ) will remain a polynomial of the same order. Thus analyticity properties of 

<72 and <73 in the extended complex plane |£| > 1 guarantee the existence of polynomial- type 

exact solutions. 

Now, consider more general initial conditions of the form 

z((,0) = E^(C,0)(C-0(0))^ +M(C,0), (51) 

J = 1 

with 7j either > 0, or <0 but not an integer. Then, it is clear that for later times, we will 
have 

= it + M(C,t ), (52) 

i 

provided M((,t ) and Ej((,t ) are required to satisfy 

Mt = ?i + (I 3 M + q 2 (53) 

E jt = 91 E h + q 3 Ej + ^ 0] > (54) 

while the location of the singularities (j(t) is determined from 

6 = “?i (Cj(t),t)- (55) 

Using the projection operator V\- on (53), it can be shown with an argument similar to that 
following (48) that if M((,t) is initially analytic in |(| > 1 with a N-th order “pole” at oo, 
then it will remain so for later times. Similar property can also be shown to hold for each 
The analyticity properties of Ej ( (,. / ) and Af(£, f) thus imply that the form of an 
initial singularity of z((,t) in |£| > 1 (including that at oo) is preserved in time. 

Moreover, from Plemelj formula applied to (32), it follows that Re I((,t ) = i n the 

limit |C| > 1 + . And since Re defines a harmonic function everywhere in |£| > 1, the 

maximum principle for harmonic functions implies that Re I((,t ) < 0 in |(| > 1. Using 
this fact and (42) into (55) it follows that 

Re Y = ~Re I(Q,t) > 0. (56) 

.WJ 

This corresponds to the statement that the singularities (j move outwards away from the 
unit circle. This property in fact transcends the restriction implicit in the decomposition 
(52); it actually holds for all singularities. To see this, note that since <7 j, q 2 and q 3 are 
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known to be analytic a priori it follows from the general theory of first order linear partial 
differential equations with analytic coefficients that any singularity of z{(, t ) must propagate 
according to (55). The fact that singularities (if present) move away from the unit circle 
implies that for a bubble with a smooth initial boundary, no finite- angled corners can form, 
the only possibility that remains open is for a zero of z c to impinge on |<| = 1 in finite time 
causing a zero-angled cusp on the interface. The solutions discussed in subsequent sections 
suggest that this is not possible when surface tension is present. However, we could not 
completely rule out this possibility for more general initial conditions. 

Now suppose that all the q/s in (51) are integers, so that z(C, 0) has only a finite collection 
of poles in j(| > 1 in addition to the simple pole at C = 0, ie > the conformal mapping is 
initially a rational function. Then it is clear from the discussion above that as long as 
the solution exists, z((,t) will remain a rational function analytic in |(| < 1. In other 
words, the problem admits exact solutions where h((,t) in (50) is a rational function. [We 
remark parenthetically, however, that closed form solutions for arbitrary nonintegral 7 j are 
not possible, because z( <,i) is not analytic at C = 1/Ci(<) inside the unit circle on secondary 

Riemann sheets.] 

The theoretical discussion in this section thus shows the close connection between the 
existence of exact solutions and the analyticity of q\, <?2 and q 3 outside the unit circle. It 
is interesting to note that the equations for the interfacial displacement in a Hele-Shaw cell 
can also be cast in a form similar to (40) when surface tension is neglected (Tanveer 1993). 
In this case, the equations also admit exact solutions in terms of polynomials or rational 
functions (see, e.g., Howison 1992 for a review). The difference, however, is that in the 
Hele-Shaw problem, when a less viscous fluid is injected into a more viscous fluid, all initial 
singularities in |(| > 1 (except that at infinity) travel towards |(| = 1 (Tanveer 1993). On 
the other hand, when the more viscous fluid displaces the less viscous one, all singularities 
move outward. In our case, the outward motion of singularities in the ( plane (away from 
the physical interface) is not affected by suction (m < 0) or injection of bubble fluid (m > 0) 
and is found to be the same for a drop (Tanveer and Vasconcelos 1994b) as for a bubble. In 
this sense, the dynamics here is very different from that of a Hele-Shaw cell. 

4. Polynomial Exact Solutions: General Results 

In this section we present a general class of exact solutions for the case in which the 
function h((, 0) is a polynomial of degree N. Such initial conditions form a dense set on the 
class of smooth initial shapes, so that any given smooth initial condition can be approximated 
by such fc.(C, 0) to any desired accuracy. The arguments of the previous section guarantee 
that h{(,t) will remain a polynomial of the same order N. Because of the outward motion 
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of complex singularities for more general initial conditions, it may be expected that the 
polynomial class of solutions t ) will be dense in the class of all analytic solutions. Hence if 
we start with a polynomial approximation /?(<^, 0) to the true initial shape, the corresponding 
polynomial is expected to remain a good approximation to the actual interface for 

later times. Thus, the solutions presented below are expected to describe the evolution for 
essentially any smooth initial bubble shape. 

Accordingly, we seek solutions of the form: 

2 (co = ^+x>,(<)e, (57) 

where the bj s are complex coefficients [and recall a(t ) < 0]. The problem now consists in 
finding a set of evolution equations for the coefficients a(t ) and bj{t ). To this end, we first 
multiply both sides of (37) by to obtain: 

2*c(c,o<?(c,o = ^(c,o{^(r 1 ,o-r 1 [/(c,o+^]^(r 1 ,o} 

+ - C^c(C,*)«c(C,0 

— [*c(C>0 + ( 2 cc(C>*)] + f/i]| (58) 

According to (18) and (25), as ( approaches zero the singular behavior of the left hand side 
of (58) is given by 


2* { (C, <)<?«,<) ~ - W - MW - ^ + 0(C). as C-0. 


C 3 • C 2 

Next consider the Taylor series expansion of I((,t ); 


(59) 


r«.0 = ^o + Et«C J . 


i - 1 


where the coefficients 7 0 and 7, are given by [see (32)] 


(60) 


/o 


Ik 


- 1 
4tt Jo 

- / 
2it Jo 


dv 


Me‘V)r 


~—iku 


k(e'V)l 


dv , k > 1. 


(61) 

(62) 


To simplify the notation we define ! 0 = I 0 + iK. Now let R((, t ) denote the right hand side 
of (58). Inserting (57) and (60) into (58), and after performing a tedious but straightforward 
algebra [made easier with the help of a symbolic-computation software such as Maple], we 
find that the singular behavior of R((, t ) as ( approaches zero is given by 

N+ 2 

fl(C,f)= ZnWr + 0(0, as (-*0, (63) 

;= i 
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with the coefficients r* of the form 


N+2-k a 

rk = -X k -(k- 1) J2 h x *+ji l < k < N + 2, 

j=o 

where the quantities Xk are given by 


Xn+2 

Xn+i 


X k 


X 2 


Xi 


abpf 

N+l-k 

abU ~ Y 3<k<N 

j- 1 

N - 1 

” jbjbj+i 
i=l 

« 2 -D N*- 

j=i 


(64) 


(65) 

( 66 ) 

(67) 

( 68 ) 

(69) 


Here star stands for complex conjugation. 

Comparing (63) with (59) and matching the terms corresponding to the double pole at 
£ = 0, one obtains the function C(t) in terms of the coefficients a(t ) and bk(t)- 



N 

X2 + X/ • 

3=0 


(70) 


Similarly, by matching the remaining singular terms we obtain the following system of ordi 
nary differential equations (ODE’s) : 


X = ™ (71) 

7T 

X t = -(k- 1) Y + «> - 3, 3<k<N + 2. (72) 

j = 0 

Note that since the area A enclosed by the curve obtained as the image of the unit circle under 
the mapping (57) is given by A = ttX,, it follows that (71) simply recovers the condition 
A = in. Equations (71) and (72) in general give a system of 2JV + 1 ODE’s from which one 
can compute the 2N + 1 parameters of the conformal mapping for specified initial data and 
external flow. In the rest of the paper we will discuss several particular cases of interest. 


5. Bubble in a linear flow 

Throughout this section we assume that the bubble has a constant volume so that we set 
m - 0. The simplest scenario is obtained when there is no imposed external flow, in which 
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case all initial shapes relax to a circle. Here however we shall consider the more interesting 
case in which the bubble is placed in a linear flow. Although the general solutions described 
in §4 can handle a large class of initial shapes, for simplicity we shall focus here on the 
case where the bubble possesses an initial circular shape. Below we present exact solutions 
describing the subsequent bubble evolution (deformation) for two specific flows: (i) pure 
straining flow and (ii) simple shear flow. 

5.1 Pure straining flow 

Here we imagine that the bubble is placed in a pure straining flow: Uo = (Qx,— Qy), 
where Q is the rate of shear. In view of (7), we thus set a 0 = 2 Q and j3 0 = u> 0 = 0. Under 
the assumption that the initial shape is a circle of radius R, one can easily verify from the 
general solutions presented above that for later times the interface will be described by a 
conformal mapping of the form 

*«>*) = ^ + K*)C, (73) 

where b(t ) is real. From (71) and (72) one finds that the time evolution of the coefficients a 
and b will be governed by the following ODE 


^ ( a fy — 2I 0 ab + Qa 2 , 

(74) 

with a(0) = R and 6(0) = 0, together with the area condition 


a 2 - 6 2 = R 2 . 

(75) 

Here the quantity 7o reads [see (61)] 


j t t T dv 


2x Jo | a 2 _j_ 5 2 — 2abcos 

(76) 


In Fig. 1 we have integrated (74) in nondimensional units for the case Q' = RQ/t = 0.75. 
We thus see that (starting from a circle) the bubbles will evolve through a series of elliptical 
shapes towards a steady solution, whose parameters Go and bo are given by the solutions to 
the equation 

Qa - 2 1 0 b = 0, (77) 

subject to the condition (75). This steady solution was first found by Richardson (1968) 
through a direct steady-state calculation. 

5.2 Simple shear flow 

Now we consider the case of a bubble placed in a simple shear flow: u 0 = (Tj/, 0), where 
r is the shear strength. Accordingly, we set o 0 = 0 and fl 0 = -lo 0 = T. If the initial shape 
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is a circle we then have, as before, that the subsequent shapes are given by the conformal 
mapping (73), but where now b(t) is a complex coefficient. According to (72), the time 
evolution is given by the following ODE 

— (ab) = -(2I 0 + iT)ab + iTa\ ( 78 ) 


plus the area condition 


|6| 2 = R 2 . 


(79) 


Here 7o reads 


/o=f r 


dv 


(80) 


{a 2 + \b\ 2 — 2ab R cos v + 2 abi sin v } 1/?2 

where b R and 6/ are the real and imaginary parts of b, respectively. 

Now defining X = ab R and Y = aft/, one readily sees that (78) yields the following system 

of ODE’s: . 


X = -(2IoX + TY) 

Y = -(2IoY + TX) + Ta 2 . 


(81) 

f'82'1 


The steady solution in this case is thus clearly given by the solution to the equations 


2l 0 b R -\-Tbi — 0 
2I 0 b I + Tb R -Ta = 0, 


(83) 

(84) 


subject to the constraint (79). [This steady solution, in a somewhat different notation, was 
also first obtained by Richardson (1968).] In Fig. 2 we show the evolution of interface towards 
the steady solution for the case = RT/t = 1. 

6. Contracting/expanding bubbles in a quiescent flow 

In this section, we suppose that the bubble is placed in an otherwise quiescent flow, i.e., 
we set o 0 = Po - u 0 = 0. For simplicity, the expansion or contraction rate m will be taken 
to be a constant throughout this section. We also assume that the bubble is symmetrical 
with respect to the x axis so that the coefficients bf s are all real. In this case, (72) in general 
gives a system of N ODE’s, which together with the area condition (71) determine the N + 1 
coefficients of the conformal mapping. Before discussing the general case, however, we will 
first consider the simpler case in which the initial shape has either elliptical or n-th fold 

symmetry. 
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6.1. ‘Symmetrical’ bubbles 


Here we seek solutions of the form 

2 (Cj t) = ~^r' + bw(t)( N , (85) 

where b N (t ) is real and assumed for definiteness to be negative. Thus, for N = 1 the bubble 
is an ellipse, whereas for N > 1 it posses (N + l)-th fold symmetry. We also note that in 
this case I k vanishes identically for 1 < k < N + 1. Using this fact on (72), it then follows 
that if the b k s are initially zero for k = — 1, they will remain so for later times, 

thus guaranteeing that a mapping of the form (85) does indeed give an exact solution. [We 
remark parenthetically however that this is not case, for example, when the bubble is placed 
in an external flow.] The evolution equation in this case is [see (72)] 

— ( ab n) = —(N + 1 )/ og 6 ; v , ( 86 ) 

while the area condition (71) after integration reads 

A(t) = 7T Ja 2 — = [T(0) + mt \ , (87) 

where .4(0) is the initial bubble area. Here Iq can be written as 


h=~ r - 

Z7T J 0 (a 


{a 2 + N 2 bx — 2Nabw cos v) 2 

as one can easily verify by inserting (85) into (61) and performing a trivial change of variables. 

Note also that a linear stability analysis of (86) readily shows that an expanding circle 
is stable whereas a contracting one is not. To analyze the behavior of the solutions above 
in more details, it is convenient to study the motion of the critical points of the conformal 
mapping, i.e., the points at which z^(C, t ) vanishes. Denoting such points by (o*., k = 0, ..., N , 
we clearly have: ( ok = [pe i2nk J 1/(;V+1) j w h ere p = a/Nb N . Using (86) and (87), we find that 
the quantity p evolves according to the following equation 


. _ p(Np 2 — 1 ) m 

p ~ WTT r + I)/ “ + ^y. • 

We also note, for later use, that in terms of p and A(t) the coefficients a and b N read 


b% = 


" N 2 p 2 N(Np 2 — 1 ) ' 

In the case of an expanding bubble (m > 0), we immediately see from (89) that p 
increases monotonically with time, so that the zeros (o k move away from the unit circle 
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(ICI = i). Hence the bubble has a tendency to become ‘smoother’ with time and asymptotes 
an expanding circle as t -> oo. An example of such a case is presented in Fig. 3 where we 
show a sequence of interface shapes for N = 1 (elliptical bubbles) and m' = 1 [see (8)]. It 
should also be noted that the results discussed in this paragraph are true in general, that 
is, for any given initial shape [described by a polynomial fe(C,0)] the bubble will approach 
an expanding circle as t -* oo. [Note, in particular, that this is true even in the absence of 

surface tension.] 

Consider now the case of a contracting bubble, i.e., m < 0. Here, of course, solutions 
can exist at most up to the time tj = A(0)/|m|, at which the bubble-fluid (‘air’) would be 
completely removed from the liquid. We will see below that in the absence of surface tension 
(r = 0), the solutions will in general breakdown before this time, due to the formation of 
cusp-singularities on the bubble surface. For r / 0, however, the solutions will always exist 
all the way up to the final time tj. 

We first consider the zero surface tension case. Setting I 0 = 0 in (89) and integrating the 


resulting equation, one finds 

„{t) = -,i + + N^i\ 

(91) 

where f = tj — t and 

m q 

7 - 2nNa{0)b N (0) > 

(92) 


Here a(0) and b^( 0) are the prescribed initial data. Thus, in this case the zeros will hit the 
unit circle at a time t c = t f - t r{N - l)/2^, since p(t c ) = 1. For N = 1 (elliptical bubbles), 
we then see that the zeros hit the unit circle at precisely the time when the bubble area goes 
to zero. Note that since a(t f ) = b{t f ) / 0 [see (90)] the final stage of the bubble in this 
case is a slit of extension 2|a(f/)| = 2yja(0)b\(0). On the other hand, for N > 1 the zeros 
impinge on the unit circle at t = t c <tj, leading to the formation of cusps on the interface 
and hence the breakdown of the solutions before the air is completely removed. 

A non-zero surface tension is expected, on general grounds, to prevent the development 
of actual cusp by providing a ‘regularization’ mechanism. Indeed, one of the advantages of 
the solutions above is that they are simple enough to allow a detailed analytical investigation 
of the regularizing effects of surface tension. For example, recalling that I Q diverges if a zero 
of 2 C (C-0 lies on the ( unit circle (|C| = 1), we immediately see from (89) that the zeros 
must approach [Cl = 1 in the limit that the area vanishes, but they cannot hit the unit circle 
while the area remains finite. Furthermore, here it is also possible to carry out an asymptotic 
analysis of the final stages of the bubble evolution, as indicated below. [A similar study has 
been recently performed by Howison and Richardson (1994) for the case of blob of viscous 
fluid with suction, where exact solutions with an analogous structure have also been found.] 
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We begin our analysis by noting [see the Appendix] that the leading-order asymptotics 
of Iq in the limit p — > 1 + is given by 


r ~ 

0 2n 


Np 2 - ll 1/2 


In (p- 1) + 0(1). 


i J 

Inserting this into (89) and writing A(t) = —mi, one finds that as p — * 1 


(93) 


n rv/ 

fj ~ — 


Np 2 -f 1 


[ t(N + 1) 

'Np 2 - 1' 

[ 2tt 

—Nmt 


1 1/2 


ln(p- 1) + 


})■ 


(94) 


The asymptotic behavior of p(t) for t -> 0 can now be calculated by balancing the most 
singular terms in (94). Here there are two distinct cases to consider: (i) N = 1 and (ii) 
N > 1. 

For N = 1 (elliptical bubbles), the contribution arising from surface tension effects [i.e., 
the first term in the right hand side of (94)] is small compared to the other terms, so that 
effectively the time evolution of the solutions is described by 

G p 2 - 1) i 


P rsj — 


as t 

t 


0, 


(95) 


which upon integrating yields 

P ~ 1 + St, as i — ► 0, (96) 

for some positive 6. Now according to (90), a linear behavior in p(t) as p — > 1 implies that 
a (tf) = b(tf ) ^ 0, so that elliptical bubbles will shrink to a slit even when r ^ 0. Here the 
effect of surface tension is simply to reduce the size of such a slit. In other words, the larger 
the surface tension, the smaller the extension of the final slit. 

For N > 1, as the zeros approach the unit circle, surface tension effects tend to slow down 
their motion, so that ‘narrow structures’ (near-cusps) are formed on the bubble surface. An 
example of this process is given in Fig. 4, where we show a sequence of interface shapes 
for N = 3 and m' = — 1, up to the formation of near-cusps. The ensuing ‘slow’ dynamics 
of the zeros, as the bubble continues to contract, can be estimated by balancing the two 
leading-order terms in the right hand side of (94), that is, 

H 1 / 2 


t(N+ 1 ) 


2tt 


N 


Nmt J 


ln(p - 1) ~ — j 


Solving this for p yields 


where 


-s/Vt 


6 = 


1 + e 


2tt 

't(N + 1) 


as t 


Nm 
' N - 1. 


0, 


i 1 / 2 


(97) 


(98) 


(99) 
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Note that, as a consistency check, one can now use (98) to verify that the left hand side of 
(94) is indeed small compared to the leading-order terms on the right hand side, as initially 
assumed in the argument above. 

To obtain a qualitative description of the asymptotic shapes as t — * 0, let us denote 
by R m ax and R min the maximum and minimum radial distance of a point on the interface 
(relative to the bubble center). From (85), one readily finds that R m ax = |« + 5| and 
Rmin = I® — 6|. Defining the deformation D ~ Rmin! Rmaxi then have 


N p — 1 
Np + 1' 


( 100 ) 


Since the deformation D goes to a constant in the limit p — ► 1 , it follows that during the late 
stages the bubble will shrink to a point through a succession of geometrically similar shapes. 
[An analogous phenomenon has also been observed by Howison and Richardson (1994) for 
the case of a blob of viscous fluid with suction.] Here the main effect of surface tension is 
to determine the time scale for when the ‘near-cusps’ first appear, that is, the greater the 
surface tension, the later the near-cusp will develop. 


6.2. ‘ Asymmetrical ’ bubbles 

We have seen above that in the case that the initial shape can be described by a conformal 
mapping of the type given in (57), surface tension effects guarantee that the solutions will 
always exist up to t = */, by which time the air will have been completely removed from the 
liquid. The situation changes considerably, however, when the bubble initially has neither 
elliptical nor n-th fold symmetry, that is, when not only b N ( 0) ± 0 (f° r N > 1) but also 
5 fc ( o) ^ 0 for some k < N. In such cases, the bubble will invariably develop a ‘thin neck 
whose width will go to zero at a time t = h < tj, after which the solutions cease to make 
physical sense. We thus refer to this process as bubble ‘breakup . In Fig. 5 we show a 
sequence of interfaces shapes leading to the breakup of the contracting bubble into three 
smaller ones. In the case shown in this figure, the coefficients b u b 3 and 6 5 were all initially 
nonzero and m! = -1. We also found a similar breakup when we take m to be proportional to 
the bubble perimeter instead of a constant, as appropriate for modeling a dissolving bubble 


gas. 

Since our solutions breakdown at the time when the two sides of the interface touch 
each other, we are unable to follow the dynamics of the ‘newborn’ bubbles. We emphasize, 
however, that no physical quantity blows up as the bubble approaches breakup. Here the 
breakdown of the solutions is caused simply by the loss of univalence of the conformal 
mapping for t > tf, (Tanveer & Vasconcelos 1994a). 


n i 
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7. Conclusions 


We have presented exact solutions for a bubble with essentially arbitrary initial shape 
evolving in a two-dimensional Stokes flow in the presence of surface tension. Our solutions, 
which are given in terms of a polynomial- type mapping function, include bubbles in a linear 
flow as well as an expanding or contracting bubble in an otherwise quiescent flow. It has 
been noted that the expanding bubble approaches a growing circle for later times, while a 
contracting circular bubble is unstable to disturbances and can lead to the formation of near- 
cusps or cause breakup before all the bubble-fluid is removed. The mathematical structure 
underlying the existence of a broad class of exact solutions has also been discussed in detail. 

As is well known, an inviscid bubble does not breakup under the action of a shear flow. 
Our results show however that breakup does occur (at least in plane flows) when another 
driving mechanism, e.g., suction of bubble-fluid, is present. While this general tendency to 
breakup is likely to be present in 3D axisymmetric flow as well, it not clear at this point 
that different parts of the interface will actually touch or merely tend to each other. Also, 
in our solutions the flow inside the bubble has been neglected since the viscosity of the inner 
fluid has been assumed to be negligible. Nonetheless, this flow is likely to be important near 
bubble breakup. The effect of a nonzero viscosity ratio thus needs to be examined in the 
future. 


One of the authors (GLV) would like to acknowledge financial support from The Ohio State University 
Postdoctoral Fellowship. This research was supported in part by DOE contract DE-FG02-92ER14270. 
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Appendix 

The asymptotic behavior of the quantity 7 0 given in (88) in the limit p -*• 1 can be 
lputed in the following way. First we use (90) to rewrite (88) in terms of p: 

| 1/2 

/oM, < 101 ) 


comr 


r \Np 2 -1 
to — 7T 


2tt [ NA(t) \ 


where 


Up) = f 


de 


{1 + p 2 - 2pcos(#)} ,/2 
We now separate this last integral into two contributions: 


f « W-jC 


de 


{(p- l) 2 + 2[1 -cos(#)]} 1/2 h {(/?- l) 2 + 4psin 2 (§)} ' 


i: 


de 


( 102 ) 


(103) 


where p — 1 << £ << 1. Here we have rewritten the integrand in both integrals above 
in a form suitable for the next step in the calculation. Accordingly, in the first integrals 
the leading-order asymptotics can be found by approximating cos# by 1 — -^0 and then 
performing a change of variable #=(/> — I) 17 - Similarly, in the second integral the leading- 
order asymptotics is obtained by simply setting p = 1. Within these approximations, one 
then has 

r e dv dv 

Io(p) ~ 

Jo 


— In 


{(p-iy + u 

e 


1/2 


+ 


I* dv 
Jt 2sin^ 


y>-i + M‘ ' 

In the leading-order the e dependence drops out from the above, thus giving 

fop) ~ ~ln(p - 1) + 0(1). 

Finally, inserting this into Eq. (101) yields formula (93). 


1 + 


1 


— In tan — e. 
4 


(104) 

(105) 


(106) 
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FIGURE CAPTIONS 


Fig. 1. Circular bubble placed in a pure straining flow with Q' = 0.75. 

Fig. 2. Circular bubble placed in a simple shear flow with P' = 1. 

Fig. 3. Expanding bubble with m' = 1. 

Fig. 4. The evolution of a four-fold symmetric bubble for m' = — 1. Note the formation of 
‘near-cusps’ on the innermost interface. 

Fig. 5. The evolution of an ‘asymmetric’ bubble leading to bubble breakup for m' = -1. 
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